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It is well known that the Poiseuille flow of a visco-elastic polymer fluid between plates or through 
a tube is linearly stable in the zero Reynolds number limit, although the stability is weak for large 
Weissenberg numbers. In this paper we argue that recent experimental and theoretical work on the 
instability of visco-elastic fluids in Taylor-Couette cells and numerical work on channel flows suggest 
a scenario in which Poiseuille flow of visco-elastic polymer fluids exhibits a nonlinear "subcritical" 
instability due to normal stress effects, with a threshold which decreases for increasing Weissenberg 
number. This proposal is confirmed by an explicit weakly nonlinear stability analysis for Poiseuille 
flow of an UCM fluid. Our analysis yields explicit predictions for the critical amplitude of velocity 
perturbations beyond which the flow is nonlinearly unstable, and for the wavelength of the mode 
whose critical amplitude is smallest. The nonlinear instability sets in quite abruptly at Weissenberg 
numbers around 4 in the planar case and about 5.2 in the cylindrical case, so that for Weissenberg 
numbers somewhat larger than these values perturbations of the order of a few percent in the wall 
shear stress suffice to make the flow unstable. We have suggested elsewhere that this nonlinear 
instability could be an important intrinsic route to melt fracture and that preliminary experiments 
are both qualitatively and quantitatively in good agreement with these predictions. 

PACS numbers: 



I. INTRODUCTION 



General motivation 



In this paper, we reconsider the classical topic of 
the stability of visco-elastic Poiseuille flow in the zero 
Reynolds number limit. From a weakly nonlinear expan- 
sion, we find that this flow is nonlinearly unstable for 
high enough flow rates. 

The first linear stability analysis of the flow of a so- 
called Oldroyd-B fluid — one of the simplest continuum 
models for a visco-elastic polymeric fluid with nonzero 
normal stress differences, characterized by a single relax- 
ation time A — was already performed almost thirty 
years ago Q • Since the subsequent careful linear stability 
analysis of Ho and Denn it is generally accepted that 
Poiseuille flow of an Oldroyd-B fluid is linearly stable, 
even though the stability is weak for large values of the 
Weissenberg number, the dimensionless quantity which 
measures the strength of polymer relaxation effects. The 
definition used in this paper for the case of cylindrical 
coordinates relevant for pipe flow is 
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where the term in the numerator is the normal stress 
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difference of the flow. For the planar geometry the index 
r has to be replaced by y, with the y-axis taken normal 
to the plates. The term in the denominator of (JJJ is 
the shear stress at the wall. For an Oldroyd-B or Upper 
Convected Maxwell (UCM) fluid, the unperturbed flow 
field v unp is simply parabolic, and we get 



Wi = 2A 



dvT p 



Or 



4^maxVi?. 



(2) 



wall 



Here t> m ax is the maximum velocity of the unperturbed 
profile, A is the aforementioned relaxation time charac- 
terizing the Oldroyd-B or UCM fluid, and R is the radius 
of the pipe. For the planar case, R has to be replaced by 
d, half the spacing between the plates. 

It is well known that visco-elastic Poiseuille flow is lin- 
early stable for well-established models like the UCM, 
Oldroyd-B and the Giesekus model Q|. However, 
there are good reasons to reconsider the nonlinear insta- 
bility of this flow configuration. First of all, Atalik and 
Keunings Q have recently presented strong numerical 
evidence tht visco-elastic Poiseuille flow is indeed non- 
linearly unstable for the UCM, Oldroyd-B and Giesekus 
model: when they injected the laminar flow field with 
a perturbation of sufficiently large amplitude, instead 
of dying out (like one would expect for a strictly sta- 
ble flow) the perturbation grew and saturated at a fi- 
nite value, resulting in a finite amplitude oscillatory flow. 
Such behavior, in which the flow is always linearly sta- 
ble but nonlinearly unstable, is also called a subcritical 
instability. A recent experiment [j| also supports this 
nonlinear instability scenario. Secondly, as we will dis- 
cuss in detail in section II Bl below, there are actually a 
lot of indirect indications from well-established results on 
visco-elastic Taylor-Couette flow that planar Couette or 
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Poiseuille flow between plates or in a tube might have 
a subcritical instability. Last but not least is the fol- 
lowing observation: if visco-elastic Poiseuille flow is un- 
stable, the nonlinear spatially and temporally oscillatory 
flow pattern it will give rise to will inevitably result in 
distortions of the flow after exiting the pipe or slit. In 
other words, if visco-elastic Poiseuille flow is unstable for 
large enough Wi, it automatically leads to an intrinsic 
route to "melt-fracture" type phenomena, a generic name 
for the fact that a polymer extrudate normally develops 
strong undulations or irregularities beyond some critical 
flow rate [E 0> E| • It is well-established that there are 
various mechanisms that can lead to such type of behav- 
ior (e.g., stick-slip behavior or exit instabilities leading to 
sharkskin type patterns 00))- Most of these can be (par- 
tially) suppressed by the proper choice of material or of 
the extruder shape and coating. Our main conclusion (in 
line with £|) that polymer Poiseuille flow is nonlinear ly 
unstable for large enough Wi implies that melt fracture 
type behavior is an unavoidable consequence of a nonlin- 
ear instability in the extruder. In this paper, we focus on 
the subcritical instability of visco-elastic Poiseuille flow 
itself. Our arguments and experimental support for this 
scenario are summarized in two recent letters 0, ITol ] , but 
the issue is complex and deserves further experimental 
and theoretical study. We will come back to it in the 
future. 



As stated, the calculations that we report in this pa- 
per give the details of the explicit nonlinear amplitude 
expansion for the stability of visco-elastic Poiseuille flow 
in the case of fixed average pressure gradient (we con- 
sider periodic modulations of the pressure, so there are 
periodic pressure modulations but the average pressure 
drop remains unchanged). We show explicitly that for 
large values of Wi there is indeed subcritical behavior, in 
other words that the flow exhibits a weakly nonlinear in- 
stability. Moreover, the subcritical behavior sets in quite 
abruptly around a value Wi c of order 5 for flow in a tube. 
This value is consistent with the value where melt frac- 
ture type behavior is normally reported to set in 
but somewhat larger than the critical value of 0.1 found 
numerically for for the related Oldroyd-B model with vis- 
cosity contrast 10~ 3 and a Reynolds number of 0.1 Q. 
We will discuss the possible origins of these discrepancies 
at the end of this paper. 



Since the nonlinear amplitude expansion that we will 
employ has been used mostly in other fields of physics, 
and since it involves some unusual subtleties, we present 
our results in some detail. For the benefit of the reader 
not interested in the derivation, we summarize our main 
results in section ITT 7 ! of this introduction. Before doing 
so, we first discuss the relation of this work with the 
Taylor-Couette problem. 



B. A nonlinear instability scenario motivated by 
the visco-elastic Taylor-Couette problem 

Just over a decade ago, Larson et al. 0] investigated 
the stability of the flow of a visco-elastic polymer solution 
in a Taylor-Couette cell, a cell consisting of two concen- 
tric rotating cylinders. They found that in this system 
the flow does exhibit a well-defined linear instability for 
some value of the Deborah number De, which is analo- 
gous to the Weissenberg number. Their calculations were 
done for the Oldroyd-B polymer model. In the same pa- 
per the predictions were confirmed experimentally by a 
series of measurements on a polymer solution which is 
well described by this model. This work was later ex- 
tended by Joo and Shaqfeh 0, , who considered the 
more general case of flow in a curved channel. In the 
limit of boundary driven flow this reduces to the Taylor- 
Couette case, while in the limit of static curved walls 
the flow is driven by a pressure gradient (so-called Dean 
flow). In all these cases, as well as in the experimen- 
tally relevant cone and plate geometry 0, 0] , the flow 
is linearly unstable at large enough flow velocities. The 
main conclusion from this line of research has therefore 
been that a linear instability occurs in visco-elastic flu- 
ids due to "hoop stresses" if the stream lines are curved 
[T7L ITsf . This is confirmed by the observation that the 
various stability calculations show that the instability 
threshold goes to infinity if the curvature of the walls 
goes to zero. For the Taylor-Couette system this is the 
limit in which the radius of the cylinders goes to infin- 
ity; for the Dean flow problem investigated by Joo and 
Shaqfeh Q this result is consistent with the weak 
stability of Poiseuille flow between two parallel planes. 

In the last few years, Groisman and Steinberg [T^.l20| 
have experimentally investigated the visco-elastic insta- 
bility in the Taylor-Couette system in detail. They find 
very good agreement with the theory for the onset of the 
instability when the polymer flow velocity is increased; 
however, an important result of their careful study is 
that the bifurcation at the critical Weissenberg/Deborah 
number is subcritical or "backward" pbl l2l| . 

Groisman and Steinberg [2(j have also presented in- 
tuitive arguments why the instability in visco-elastic 
Taylor-Couette flow at low Reynolds numbers is subcrit- 
ical. In their arguments the fact that the flow is oc- 
curring between curved walls is playing essentially no 
role: the curvature of the walls is necessary to make the 
flow linearly unstable, but the nonlinear positive feed- 
back mechanism that in their picture gives rise to the 
subcritical behavior does not rely on the curvature. This 
suggests the scenario sketched in Fig. ^ The stable non- 
linear flow behavior, indicated in the figure by the full 
lines, is quite independent of the precise radius or cur- 
vature of the cylinders, assuming the "gap" (distance) 
between the cylinders to be constant. The points where 
the curves touch the horizontal axis correspond to the 
linear instability threshold. These points therefore are 
strongly dependent on the curvature of the cylinders: as 
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Weissenberg number 

FIG. 1: Qualitative sketch of the subcritical behavior of the 
visco-elastic instability in low Reynolds number polymer flow 
in a Taylor-Couette cell, and of our proposal that the bifur- 
cation curve for the case of Poiseuille flow should be thought 
of as the limit in which the radius in the Taylor-Couette cell 
goes to infinity. See text. 

indicated in the figure, the results of the theoretical anal- 
ysis 0, Q, 0| imply that they shift to the right as the 
radius R of the inner cylinder increases, and that there 
is no linear instability in the limit R — » oo. 

The continuity found by Joo and Shaqfeh [T3. Il4| 
in going from the Taylor-Couette flow to the case of 
Poiseuille-like Dean flow between fixed curved planes, as 
well as the intuitive arguments of Groisman and Stein- 
berg [20] strongly suggest that the scenario for Poiscuillc- 
flow through a channel will be very close to the one 
sketched in Fig. ^ for the limit R^ — > oo: the unper- 
turbed flow is linearly stable for any Wi, but nonlinear ly 
unstable as if there is a subcritical bifurcation at Wi = oo. 
Moreover, just like the rightmost curve for R4 — > 00 ap- 
proaches the horizontal axis rapidly for sufficiently large 
Wi, implying that the threshold for the nonlinear instabil- 
ity is small at large Wi, we expect that the threshold for 
the subcritical-like nonlinear instability of visco-elastic 
Poiseuille flow is small for sufficiently large Wi. The work 
done in this paper, summarized in the next subsection, 
fully confirms this expectation. 

It is of interest to note that the scenario we propose 
here for visco-elastic fluids in the zero Reynolds num- 
ber limit has strong similarities for large Reynold num- 
ber planar Poiseuille flow of Newtonian fluids: planar 
Poiseuille flow of ordinary fluids becomes linearly unsta- 
ble at a Reynolds number of 5772, but in practice the 
flow becomes unstable at much lower Reynolds numbers 
of IH, |3] . Thus the transition is subcritical. The 
scenario that has emerged is that the nonlinear branch 
extends down to Re w 2500 for two-dimensional pertur- 
bations and to Re rj 1000 for three-dimensional flows, 
and in fact already over thirty years ago [3, |2(| |27], |28| , 
an amplitude expansion for planar Poiseuille flow in the 
spirit of ours was instrumental in showing that the insta- 
bility at Re = 5772 is subcritical. The scenario we pro- 



pose for visco-elastic flow is even closer to the one that is 
found for the transition to turbulence in Poiseuille flow 
of Newtonian fluids in a pipe: although the flow is also 
linearly stable for any Re, the flow is nonlinearly unstable 
for Re > 1000 with a threshold which decreases as Re -7 
for Re — ► 00 [2|j, very much like we suggest in Fig. 
with the dashed line labeled R4 — > 00. The transition 
to turbulence of Newtonian fluids therefore shows that 
whether or not there is a true linear instability is not 
very relevant in practice, and we believe that this is true 
for visco-elastic flows too. 

Let us conclude this subsection with the following ob- 
servation. As noted in the preceding paragraphs, it 
is generally accepted that when the stream lines are 
curved, visco-elastic flow becomes linearly unstable at 
sufficiently large flow rates. In accord with this visco- 
elastic Poiseuille flow is linearly stable since its stream 
lines are straight: infinitesimal perturbations must and 
do decay. However, following this same line of reasoning 
it is also very natural for visco-elastic Poiseuille flow or 
planar Couette to be nonlinearly unstable: in the pres- 
ence of a finite perturbation, the stream lines are curved, 
and consequently we expect the flow to be unstable to 
any finite perturbation for sufficiently large flow rates. 
This is the essence of the subcritical instability scenario 
sketched above. 



C. Summary of our main result for Poiseuille flow 
in the UCM model 

In this paper we will consider the limit that the 
Reynolds number Re 

is negligible. Here R is a characteristic length scale (d for 
the planar case, the radius R for the case of a cylindrical 
tube), v a characteristic velocity, p the density and 77 the 
viscosity of the fluid. As Re measures the strength of the 
inertial terms with respect to the viscous terms, in the 
limit Re J, we can ignore the nonlinear convective terms 
in the momentum (Navier-Stokes) equation, and make 
a quasistationary approximation in which the temporal 
derivatives in the momentum equation are neglected. 

As stated before, as the constitutive equation for the 
polymer fluid we take the so-called UCM (Upper Con- 
verted Maxwell) model |l|, which expresses the stress 
tensor t of the polymer fluid in terms of the shear tensor 
through 

?+ Ar (1) = -ry(Vu + (W) 1 ), (4) 

where "the upper convected derivative" fvx) is given ex- 
plicitly in Eq. (|19fl below, and where A is the parameter 
with the dimension of time that characterizes the UCM 
model. 
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In this paper, we consider a perturbation of the veloc- 
ity and stress fields with single wavenumber k along the 
direction of the flow, and with amplitude A(t), i.e. 



perturbed fields oc A(t)e + c.c 



(5) 



where c.c. means complex conjugate. Then, in an expan- 
sion in powers of A, we determine to lowest nonlinear 
order the equation for A, 



A A 

^- = - iL o(k)A + c 3 \A\ 2 A. 



(6) 



To linear order in A this equation simply reproduces the 
term iuj(k) of the dispersion relation of a single mode 
gifcz-jwt. although we have to redo the linear stability 
analysis in order to proceed to the nonlinear term, in 
principle this term is already contained in the analysis 
by Ho and Denn Q. In particular, since we know that 
every mode k is linearly stable, Imw(fc) < for all k. 
The essence of our analysis consists of calculating the co- 
efficient C3 explicitly. In particular the real part of C3 is 
of importance: if the real part Rec3 < 0, then the non- 
linear terms increase the damping of the amplitude and 
the unperturbed state is, within this approximation, not 
only linearly, but also nonlinearly stable. On the other 
hand, if Re C3 > 0, then the nonlinear term promotes the 
growth of the amplitude, and in particular amplitudes 
satisfying 



\A\ >A C = 



llmuj(k) 
Rec3 



(7) 



grow without bound. Hence, in this approximation A c 
defined above constitutes the critical amplitude of the 
perturbation beyond which the flow in nonlinearly unsta- 
ble. 

In our analysis, we do find that indeed for sufficiently 
large Weissenberg number Re C3 > 0; we then determine 
the value of k for which A c is smallest, and take this as 
the critical amplitude for the nonlinear flow instability. 
The value of A c obtained this way from our analysis is 
plotted as a function of Wi in Fig. [21 for the planar case 
and in Fig.|3for the cylindrical case. Our normalization 
is such that A is the ratio of the maximum perturbation 
in the shear rate at the wall over one wavelength, divided 
by the unperturbed shear rate , 



max [d Sv z /dy] 



dv z p /dy 



(8) 



wall 



As we see from Figs. |2 and 01 that the overall behavior of 
the critical amplitude required to trigger the instability 
is in accordance with the picture as suggested in Fig. I: 
for i?4 — * 00, the planar limit, the threshold amplitude 
is expected to get increasingly small as one increases Wi, 
as indeed it does. 

In Fig. 21 we show the velocity profile Sv correspond- 
ing with the linear eigenmode with wavelength A = 1.7 R 



in the planar case. This wavelength is close to the one 
with the lowest instability threshold. The roll-type struc- 
ture of the flow-profile (which we also find in the planar 
case, see section III) has very much the same structure 
as the one in the Taylor-Couette cell that according to 
the arguments of Groisman and Steinberg [2(j underlies 
the subcritical instability in that case — this confirms 
that essentially the same mechanism is responsible for 
the subcritical instability in Poiseuille flow. 

An important feature of our results which is of prac- 
tical importance is that the nonlinear instability sets in 
quite abruptly: below some critical value Wi c of the Weis- 
senberg number, the flow in the cylindrical case is within 
our approximation also nonlinearly stable as Re C3 < 0, 
while above Wi c the critical amplitude, especially the 
critical amplitude for the shear stress at the wall, drops 
rapidly to a small value (as we shall see, the case of pla- 
nar Poiseuille flow is slightly different). We can make 
this more precise as follows. In the general scenario, the 
nonlinear branch where the flow-profile is nontrivial and 
of the form sketched in Fig. 01 ends at a so-called saddle- 
node bifurcation point — this is the point marked with 
a dot in Fig. ^ where the unstable dashed branch and 
the stable branch, indicated with a full line, meet. For 
control parameters below the value corresponding to the 
saddle-node bifurcation, the nontrivial flow pattern is not 
dynamically stable anymore. Now, as we extend our ex- 
pansion only up to the first nontrivial order in the am- 
plitude, we are unable, from our expansion, to precisely 
locate the saddle-node bifurcation point. However, it is 
reasonable to associate the approximate location with the 
point where Rec3 = 0. In this approximation, our esti- 
mate for the saddle-node bifurcation, and hence for the 
possible onset of the modulations in the flow profile and 
hence in the distortions of the extrudate, is Wi^ yl « 5.2. 
Note also that the threshold value for the nonlinear insta- 
bility drops quite sharply when Wi increases beyond this 
value: Assuming even careful experiments can not avoid 
perturbations of the order of a few percent, we would 
expect to actually see the nonlinear state at Wi values 
somewhat larger than \N\^. yl « 5.2 in the cylindrical tube. 
This is quite consistent with the Weissenberg value Wi c 
where according to the literature 0, 0, 13 extrudates typi- 
cally exhibit undulations and deformations. In the planar 
case the value of WiJ? 1 is somewhat less sharply defined, 
according to our results, but for all practical purposes 
it appears to be somewhat smaller than Wi^ yl . We will 
discuss the competition between the unperturbed lami- 
nar flow profile and the nonlinear profile further in the 
concluding section HV1 

A few remarks concerning Figs. |3 and |21 which consti- 
tute the main result of our analysis, are in order: 

(i) Of course, our analysis is only based on an expan- 
sion to lowest nontrivial order in A] one may therefore 
wonder how much higher order nonlinear terms affect 
the result for A c . Clearly, for small Wi where A c be- 
comes of order unity, our results should not be trusted 
quantitatively: higher order (quintic) terms will change 
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FIG. 2: The critical amplitude and the critical stress for the 
case of planar Poiseuille flow of an UCM fluid, as determined 
from the weakly nonlinear expansion in this paper. 
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FIG. 3: The critical amplitude and the critical stress for the 
cylindrical tube. Note that the critical values of T2s/t^ p 
have been multiplied by 10 so as to be able to use the same 
scale. 



the answer (for the planare case they may even give 
rise to the existence of a true saddle-node bifurcation 
point). However, for large enough Wi, our results can be 
trusted quantitatively. This is because the linear damp- 
ing term Imw becomes arbitrarily small at large Wi (it 
varies roughly as Imw ~ l/W\), so that then the ampli- 
tude expansion for A c becomes nicely ordered. 

(ii) Since we have only determined the cubic non- 



FIG. 4: Plot of the velocity field Sv corresponding to the 
linear eigenmode of the cylindrical geometry with wavelength 
A = 1.7 R at Wi = 8. The basic flow profile is in the horizontal 
direction. 



linearity in the equation, we can not say anything about 
the finite amplitude at which the instability will saturate. 
It is in fact questionable whether the saturation ampli- 
tude can be determined from any perturbative method. 
One hint that this might be possible to a reasonable ap- 
proximation comes from the experiments of Bertola et al 
[lfl Y\A\ , which indicate that the amplitude of the pertur- 
bations of the extrudate just above Wi c is rather small. 

(Hi) The scale on the vertical axis is not arbitrary. In 
Figs. EJand |21 we have plotted the size of the shear rate 
perturbation normalized to the shear rate at the wall in 
the unperturbed case. Using the equation 



max[()7 



unp 



- \C{1)A C 



(9) 



wall 



which holds both for the cylindrical case and the planar 
case(with r replaced by y), with C(l) a numerical con- 
stant defined in Appendix [Bj we show in Figs. [21 and |31 
the ratio of the perturbed shear stress at the wall over 
the unperturbed shear stress at the wall, beyond which 
the flow is unstable. Note also the steep drop of the curve 
for Wi just above Wi c : for all practical purposes the tran- 
sition is quite sharp. 

(iv) Our analysis also yields an idea of the value of 
k of the mode with the smallest critical amplitude. For 
large enough Wi that our analysis can be trusted, we find 
typical values about twice the diameter of the slit or the 
tube both in the planar case and in the cylindrical case. 
A precise comparison has to be based on analyzing the 
frequency of the flow distortions measured at a fixed po- 
sition, however, since the flow distorts upon exiting the 
die — see section llVl 



D. Outline of the paper 

We present our nonlinear analysis in some detail for 
various reasons. First of all, an amplitude analysis is nor- 
mally used for problems with a true linear instability; the 
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use of the method for cases like this one without a true in- 
stability involves some unusual nontrivial subtleties. Sec- 
ond, such an expansion makes use of left eigenvectors of 
the linear operator, whose behavior and boundary con- 
ditions are quite intricate and worth discussing. Third, 
the analysis we introduce may be of relevant to other 
visco-elastic flow problems as well. 

In order to perform our weakly nonlinear analysis, we 
first have to reanalyze the linear stability problem. The 
essential results are summarized in section [H] In order 
to facilitate comparison with the earlier work by Denn 
and coworkers @ Q , we write the equation for the sta- 
bility eigenmodes in terms of the stream function, which 
satisfies a fourth order linear differential equation. The 
coefficients that we have (re) derived for this differential 
equation are given in appendices, and are the same as 
those given in the appendix of 2]. A feature not dis- 
cussed in the earlier work, however, is that there are var- 
ious eigenmodes with different symmetries in the verti- 
cal direction for the planar case. In section ITTT1 we first 
rewrite the linear eigenvalue problem in a form which is 
closer to the one usually found in derivations of an am- 
plitude expansion [3(], |^|. This is then followed by a 
discussion of the derivation of the cubic nonlinearity in 
the amplitude equation. A somewhat special feature of 
our approach that should be kept in mind is that nor- 
mally amplitude expansions are used for an expansion 
around a true bifurcation point where a particular mode 
loses stability. Here the relevant linear modes are always 
weakly damped. This gives rise to some slight differences 
in the formulation. 



in the planar case and 

He = 0, , v r>z = v r , z (r,z). 



(11) 



in the cylindrical case. The first step is the derivation of 
an expression for the unperturbed flow: 

« U "P = (0, o, vT p (y)), « unp = (0, 0, < np (r)). (12) 

The next step is the perturbation of the unperturbed 
flow: 

v = ^ + (0,8v y (y),Sv z (y)y {kz - UJt \ (13) 
v = ^ n P + (Sv r (r),0,Sv z (r)), (14) 

m = T%» + 5T ij e$ k '- ut K (15) 

We will keep the terms linear in STij and 5v and write 
the Srij and 6v z in terms of dv y (planar) or 5v r (cylindri- 
cal) and its derivatives. This will give us a fourth order 
equations for Sv y or Sv r . with boundary conditions. Solv- 
ing this equation yields the dispersion relation u{k); the 
imaginary part of u> determines the stability of the flow. 



B. Equations and boundary conditions 

1. The equations 

The flow is taken to be incompressible, so that the con- 
servation of mass equation becomes the incompressibility 
condition 



V -v = 0. 



II. LINEAR STABILITY ANALYSIS 

After having formulated the problem for the planar 
case in the first subsection, we will summarize the main 
steps of the linear stability analysis in the second sub- 
section. Detailed expressions for the coefficients are rel- 
egated to to appendices. The numerical results for the 
dispersion relation of the linear modes are presented in 
the last subsection. 



A. Formulation of the problem 

We would like to investigate the linear stability of poly- 
meric flow between two plates, separate by a distance of 
2d, and through a cylinder of radius R. The direction 
orthogonal to the planes will be our y-direction, the di- 
rection of the unperturbed flow will be the z-direction 
- the advantage of this convention is that both in the 
planar and in the cylindrical geometry the mean flow is 
in the z-direction. From now on we take the flow two- 
dimensional by putting 



0. 



(10) 



(16) 



Moreover, as was already explained in the introduction, 
since we are interested in the zero Reynolds number limit, 
the Navier-Stokes equation, expressing conservation of 
momentum, reduces to the linear equation 



-Vp- V • f = 0. 



(17) 



We will take the curl of equation (|17fl to eliminate the 
pressure. This leaves us with an equation for the com- 
ponents of the stress tensor. The UCM model describes 
the stresses in the polymeric fluid [1|: 



Af (1 



(i) 



-ry(Vu+ (W) f ), 



(18) 



where the stress tensor f^) is defined in the following way 




= §^-(^) t -^-(^)> (I 9 ) 



Di = dl +V - V - 



(20) 



Eq. (|18fl illustrates that the UCM model is characterized 
by one time constant A, which models the polymer relax- 
ation time. The "upper convected derivative" rm is the 
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simplest frame-independent formulation that implements 
this Q. 

The UCM constitutive equation only models normal 
stress effects, no shear thinning. This is illustrated by 
the well-known fact that upon using the Ansatz (|12|1 for 
v unp , we find that the steady flow profile of a UCM fluid 
is still parabolic, 



1 



planar case, (21) 



where v max is the maximum velocity at the center line. 
Furthermore, for the stress tensor of the basic parabolic 
profile, one finds for the planar case simply 



T unp _ unp _ 

' yz zy 



unp 



-2r]X 



<9< np (y) 

dvT p (y) 

dy 



(22) 



(23) 



All other elements of r unp vanish. Note that r"" p is 
nonzero and proportional to the square of the shear — 
this is the normal stress effect. 

The results in the cylindrical case are: 



vT V {r) = %ax 



1 



and 



.unp 



.unp 



-2r?A 



-'/- 



<9< lp (r) 
dr 

<9< np (r) 
dr 



(24) 

(25) 
(26) 



The first step in a linear stability analysis is to linearize 
equation (T5)| to obtain expressions for the 6r which are 
linear in 5v. At this stage, it is useful to introduce the 
stream function <!>. Generally, the stream function is in- 
troduced by writing (for the planar case) 



<9$ 

dz ' 



d<P 

dy ' 



(27) 



so that the incompressibility condition l|16fl is satisfied 
automatically. For linear perturbations of the form 1|13|) 
we simply have $ = 4> e^ fcz_wt ^, so that 



5v y (y) =ik(f>(y), Sv z (y) 



d<j){y) 
dy 



(28) 



The stream function is slightly more complicated in the 
cylindrical case: 



—4>{r) = Sv r (r), 
r 



1 d(j)(r) 



5v z (r). (29) 



r dr 

It is also convenient to introduce dimcnsionless variables 



plane: z 



z/d 

V 

d' 



k = kd, 
i=tv max /d, (30) 



for the planar die, and 

Cj — 
cylinder: z 



ljR 

^max 

z/R 
r 

R' 



k = kR, 
i = tv max /R, (31) 



for the cylindrical die. For notational simplicity, we will 
rename k =>■ k, z z and t t and uj => u; it just means 
that we have to remember that all lengths are henceforth 
measured in units of d or R, inverse lengths in units of 
1/d or 1/-R, and times in unites of d/v max or R/v max . 
The imaginary part of u) determines the linear stability 
of the flow; the flow is linearly stable if Imw < 0. In 
order to facilitate comparison of our explicit expressions 
for the linear equation, given in appendix El with the 
expressions of Rothenberger et al. |2( , we also introduce 
their dimensionless number S 0, 



Wi 



(32) 



Using definition (|28|l we obtain equations for the St in 
terms of <j>{y). Equation (|16fl is always satisfied because 
of definition l|28|l , equation (|17|) will give us a fourth order 
equation for <f>(y) of the following form: 



where 



0"" + fa<j>'» + + + M = 0, (33) 



fl, =&(*:, S,w;0. (34) 



We used the symbolic manipulation program Maple to 
find the explicit expressions for the these expressions 
are the same as those given by Rothenberger et al. 
and can be found in the Appendix. 



2. Boundary and symmetry conditions, planar case 

The aim of the linear stability calculation is to deter- 
mine u>(k) for fixed k and S, such that <f> satisfies the 
usual stick boundary conditions: 



0. 



at ( = ±1. 



For <j) these boundary conditions translate into 

<H£ = ±i) = o, 0'(c = ±i) = o. 



(35) 



(36) 



Note that we have four boundary conditions, and a fourth 
order equation. At first sight, one might think that there- 
fore the equation might have unique solutions for any 
Lu(k). However, because of the linearity of the problem, 
if we find a solution <?!>(£), C<fi(t;) with C an arbitrary 
complex constant is a solution of equation (|33|l as well. 
We eliminate this arbitrary degree of freedom by setting 



(37) 
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Since this eliminates two trivial degrees of freedom which 
do not affect the solution, it is now clear that for a given 
k and S, one or more unique branches of the complex 
quantity u)(k) will be fixed by the differential equation 
for <fi. 

Because of the vertical symmetry of the problem, 
eigenfunctions will either be asymmetric or symmetric. 
In the first case the boundary conditions on the center- 
line are 



asymmetric: v y — 0, 



u;'=0 at f 



0, 



which implies for cj>(y) 

asymmetric: cj) — 0, tfi" = at £ = 0, 
and in the second case the conditions are 

symmetric: v' y = 0, v' y " — at f = 0, 
which implies 

symmetric: <f>' = 0, 4>"' = at £ = 0. 
We will explore the stability of both profiles below. 



(38) 
(39) 
(40) 
(41) 



3. Boundary conditions, cylindrical case 

We will take the usual stick boundary conditions v r — 
0, v z = at £ = ±1. In the center we require v r = 0, 
v z finite. In terms of <j> we have <j> = Q, (/>' = at £ = 1, 
and we choose </>"(l) = 1. We have to be careful in the 
origin, because of the 1 /£ terms in the equations for the 
/3j. The boundary conditions on v imply (j) = 0, 4>' = 
at £ = 0. In order to avoid numerical problems, we use a 
polynomial <£(£) = a 2 £ 2 / 2 + a 3 £ 3 /6 for £ 6 (0,0.01). At 
£ = 0.01 we match both solutions by imposing equality 
of <b and the first three derivatives. 



4- Summary of the linear stability problem 

In summary, we have to solve equation (|33[) with the 
condition l|37|l and the boundary conditions (|36() and 
(|39|l for asymmetric modes or with the boundary con- 
ditions H3ti|) and (|41|l in the case of symmetric modes. 
This means that we have to find the two complex pa- 
rameters 0"'(1), u>(k), such that the conditions on the 
centerline are satisfied. In the cylindrical case we have 
to solve equation (|33|l with the boundary conditions dis- 
cussed above. This means that we have to find the four 
complex parameters </>"'( 1), ui, a 2 and a 3 such that both 
solutions can be matched together. 



Numerical results 



The planar case 



We have used a shooting program [3^J to find the right 
parameters and to construct the eigenmodes correspond- 



-0.2 



3 -0.4 
E 



-0.6 



-0.8 















•0 000^ 








—•— Wi=4 
Wi=6 
Wi=12 
Wi=20 



















































1.2 



1.4 



1.6 



1.8 



FIG. 5: Dispersion relation for Imw = — e of the weakest 
damped asymmetric mode for Wi = 4, 6, 12 and 20 in the 
planar case. 



ing to the two classes of boundary conditions. We present 
the results below for the range of /c-values for which our 
program converged essentially to arbitrary precision. For 
larger values of k the convergence becomes poorer, but 
since these larger values do not appear to be relevant for 
the nonlinear analysis of section lnTl we content ourselves 
with reporting simply the range where sufficient precision 
could be reached. 

Fig. [S] shows the magnitude of the imaginary part of 
the eigenvalue u) of the asymmetric mode as a function of 
k and for various Wi. For all values of Wi, Imw < 0. As 
the temporal behavior of the modes is as e _et , this con- 
firms that the flow is linearly stable. Especially for large 
Wi, e is very small, which means that the flow is only 
weakly stable. This is the reason that we also introduce 
the dimensionless temporal decay rate e = -Imw, so 
that the positive quantity e is a small quantity for suffi- 
ciently large Wi. This is very important, as the smallness 
of e will allow us to use an amplitude expansion in the 
next section: this expansion is based on an adiabatic ap- 
proximation for the growth of the amplitude relative to 
the intrinsic oscillation of the waves with frequency Rew. 
Since the intrinsic frequency is of the order of unity in 
our dimensionless units, the condition for the amplitude 
expansion to work is that s <C 1. 

We found a second asymmetric mode close to the first 
one reported in Fig. [3] This mode is slightly more 
damped than the first one. The two asymmetric modes 
are shown in Fig. |SJ The symmetric mode lies between 
the two asymmetric modes, for, quite surprisingly, we 
have been able to show analytically that the symmetric 
mode obeys 



Im u>(k) 



4 

Wf 



(42) 
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FIG. 6: The dispersion relation for the two different linear 
asymmetric eigenmodes at Wi = 8 in the planar case, and 
the symmetric mode whose damping rate is given exactly by 
Eq. tHJ. 



so that for the the Weissenberg number Wi = 8 in Fig. El 
the damping rate of the symmetric mode is e = 0.5. For 
this mode the phase velocity Reu>(k)/k is extremely close 
to 1, but the precise value has to be determined numer- 
ically. Since the damping rates of all the three modes 
are very close, the above analytical result nicely shows 
that the damping rate of the modes becomes arbitrarily 
small for sufficiently large Wi = 2S. Therefore, our am- 
plitude expansion becomes self-consistent for sufficiently 
large Wi. 

Once the eigenmode has been obtained numerically, 
one has the velocity and shear stress fields as a function 
of the coordinates. Fig. [S] of the introduction illustrates 
perturbed velocity field in the planar case for Wi = 8 and 
fc=1.5. 

Our nonlinear analysis in the next section will be based 
on using the asymmetric mode which is the weakest 
damped. The reason for not basing our expansion on 
the symmetric mode is twofold: 

(i) The symmetry of the velocity v y of this mode is 
such that it corresponds to a type of undulation mode 
between the plates, which does not generalize easily to 
the case of a cylindrical tube. 

(ii) The result that the eigenvalue u>(k) = (1 + S)k — 
2i/S for this mode, with S a very small real quantity, im- 
plies that the factor c(x) which appears at various places 
in the equations for the (see Appendix |B|) . becomes 
very large near the center. Consequently, the compo- 
nents of r become very large, almost singular, near the 
center. 

(Hi) The asymmetric mode is the one least damped. 
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FIG. 7: Dispersion relation for the cylindrical case for Wi 
4, 8, 10 and 12. 




FIG. 8: Plot of the velocity field Sv corresponding to the linear 
eigenmode mode of the planar slit geometry with wavelength 
A = 4n/3 at Wi = 8. The basic flow profile is in the horizontal 
direction. 



2. The cylindrical case 

We now turn to the numerical results for the cylindri- 
cal case. Fig. [7| shows the magnitude of the imaginary 
part of the eigenvalue u> as a function of k and for various 
Wi. As in the planar case, all values of Wi, Imu < 0. As 
the temporal behavior of the modes is as e~ luJt , this again 
confirms that the flow is linearly stable; the decrease of 
e with increasing Wi again confirms that the linear sta- 
bility becomes less and less with increasing flow velocity, 
roughly as e ~ 1/Wi. 

In Fig. 01 we already showed the velocity profile Sv cor- 
responding with the linear eigenmode with wavelength 
A = 47r/3 in the cylindrical case. Fig. PJ] confirms that the 
same type of behavior is found in the planar case; the 
main difference is that in the cylindrical case, the mode 
is more confined close to the center of the tube (see also 
Fig. [5J than in the planar case. This is consistent with 
the fact that, as we will see in Section IV, the nonlin- 
early most unstable mode has a shorter wavelength in 




FIG. 9: The data of Fig. |I| plotted differently. Here the com- 
ponent 8v z is shown in a three-dimensional plot; it illustrates 
how the perturbation is largest at the center of the tube. 

the cylindrical tube than in the case of the planar slit. 
In fact, the roll- type structure of the planar flow-profile 
is more evenly distributed over the whole cross-section, 
its qualitative appearance is closer than the cylindrical 
one to the flow pattern in the Taylor-Couette cell that 
according to the arguments of Groisman and Steinberg 
20] underlies the subcritical instability in that case. 

III. NONLINEAR ANALYSIS 

We will start with an outline of the method that we 
use in the first section. As we shall see, in important 
ingredient that we need to obtain to determine certain 
solvability conditions is the solution of the linear adjoint 
operator. We will discuss the derivation of the adjoint 
problem in the sections tlll fjl and llll CI We finally discuss 
the numerical results in section ITlI Dl 

A. Outline of the structure of the amplitude 
expansion 

In order to perform the nonlinear analysis it is conve- 
nient to rewrite the equations in a vector notation. We 
shall develop the framework for the planar case, and then 
indicate the changes in the cylindrical case at the end of 
this section. 

Since we only consider perturbations which are inde- 
pendent of the coordinate x in the planar case, the stress 
components t xx , r xy , and t xz are all zero, also in the 
nonlinear regime. We can therefore capture the nonzero 
fields in a five-component vector V whose components 
are 5v y , 5v z , 5t22, ST23 and ST33. Furthermore, we rewrite 
equations (|16fl . (|17|) and i|18|) in the following form: 

CV = N(V, V), (43) 
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where C is the linear operator associated with the linear 
problem of section [Q] and N(V, V) contains all the non- 
linear terms. Thus the linear problem is simply CV = 0, 
and N(V, V) contains only terms from the left hand side 
of equation 1|18|) : for the analysis below it is convenient 
to take N symmetrized, so that N(Vi,V2) = N(V 2 ,Vi). 
Explicit expressions for N(V, V) and C can be found in 
Appendix ^ 

Normally, an amplitude expansion is performed about 
the critical point where one of the modes is marginal 
(i.e., neither grows nor decays) for a particular value of 
the wavenumber k [3fjl l8l| . This critical wavenumber 
corresponds to a maximum of the linear dispersion rela- 
tion. In the present case, the situation is not quite like 
this: the dispersion relation does not have a clear maxi- 
mum, and the linear modes are always weakly damped. 
This difference is in practice not a great problem. First 
of all, we will not be interested in spatial variations of 
the envelope; instead, we just pick a particular mode k, 
and determine the important one for our analysis at a 
later stage (through the requirement that the threshold 
amplitude be minimal). It is therefore not necessary to 
expand about a maximum of the dispersion relation. Sec- 
ondly, the amplitude expansion is based on an adiabatic 
decoupling of the fast and slow scales. In this problem, 
the fast scale is the period of the modes and the slow 
time scale is associated with the linear decay time of the 
modes. We saw in section ITT1 that in dimensionless units, 
the frequency of the modes is of order unity, while the 
damping rate e becomes much smaller than 1 for large 
enough Wi. Thus, there is indeed a separation of times 
scales, and this is essentially all that is needed for the 
analysis below. 

In fact, for readers experienced with amplitude expan- 
sions, we could essentially go ahead pretending we are 
expanding about a true critical point where the growth 
rate of one of the modes is zero, and then at the end add 
the linear damping term by hand. Nevertheless, we pre- 
fer to formulate the analysis more carefully by keeping 
the damping term as it stands. 

A linear eigenmodc of the equations is of the form 

V {£,z,t;k,w) = V (Z;k,u;)e^ kz -^ 

= V (£]k,w) eC***-*"'*) e- £t , (44) 

where we have introduced u> r = Hecu. Of course, all 
components of Vq can just be obtained directly from the 
results of the linear stability analysis of section [n] 

In an amplitude equation formalism, one can in general 
also allow for spatial variations of the amplitude on a slow 
scale. As we already remarked above, here we confine 
the analysis to the temporal evolution of a single mode 
k and its harmonics. We do so for two reasons: First of 
all, it simplifies the analysis, secondly, our main goal is to 
determine whether there is a weakly nonlinear instability. 
Anticipating that we will find that there is one, there is 
then not much to be gained in allowing for slow spatial 
variations. 
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Our aim thus is to study the weakly nonlinear evolu- 
tion in time of the amplitude A of a single mode of the 
form 

A(t) V (& k, u) eC**-***) + c.c. (45) 

Of course, such a single mode is not really a solution of 
the (weak) nonlinear equation, but the corrections are 
automatically accounted for in the amplitude equation 
formalism and are slaved to 1(450. We also note that be- 
cause of the normalization 1(37(1 . <fi"(l) = 1, we have 

max[5^(0|waii] = 2|A| (46) 

which, because in dimensionless units the shear rate at 
the wall in the unperturbed flow field equals 2, leads im- 
mediately to the the relation JSJ between the maximum 
relative shear rate distortion at the wall and the ampli- 
tude A. Likewise the relation © between the amplitude 
A and the maximum shear stress distortion at the wall 
follows from the explicit linearized dynamical equations. 
These identifications are important in translating our re- 
sults to real values, as was already noted in section I.C 
and Fig. H and 

We want to know in particular whether the amplitude 
A will grow or decay in time when the nonlinearities are 
taken into account. In the amplitude equation formalism, 
A(t) is expected to depend only on the slow timescale at, 
but we have made this explicit, for reasons that become 
clear below. We have to keep track of the fact that the 
true zero mode of our linear operator C is Vo as given 
in 144|) . This term includes the weakly damped expo- 
nential factor e _£t , and this term is not explicit in 1)45(1. 
The amplitude equation now proceeds by constructing 
the weakly nonlinear solution by writing |30l l31| 

V = eW s +eV 1 +eiv 2 + --- , (47) 

where V S is the real quantity defined as the sum 

V S =B(T)V + B*(T)V *, (48) 

and where in writing B(T) we have anticipated that the 
amplitude B varies on the slow timescale T = et. Note 
that the real combination V S is needed because the vector 
V refers to real physical fields. Once we have derived 
the equation for B, comparison with (|45|) shows that we 
should make the association 

A(t) e?B(T)e- et . (49) 

The other terms in 1(47(1 are then precisely the corrections 
to the dominant mode which are necessary to construct 
a weakly nonlinear solution. 

In the subsequent analysis, we have to keep in mind 
that Vo, with the temporal factor included, is the true 
zero eigenmode of the linear operator C: The temporal 
derivatives in the linear operator in fact work explicitly 
on both temporal exponential factors in the expression 



((44(1 for Vq . When substituting the form ((47(1 in the equa- 
tions, we then also have to account for the time deriva- 
tives of B(T). For a product term of the form B(T)e~ lu ' t 
we then have 

d(B{T)e- lut ) , Be- luJt „ , t dB(T) , % 

which we can simply summarize by making the substitu- 
tion dt — > dt+ sdx in all the time derivatives in the linear 
operator. Since only first order time derivatives enter the 
linear operator C, we get an expansion of C of the form 
C = C + z£t- 

The amplitude expansion now proceeds by substituting 
the expansions for C and 1(47(1 for V in the nonlinear 
equation 1(43(1 and collecting the terms order by order in 
e. For the first three orders in e we get: 

order : C V Q S = 0, (51) 

order e 1 : CqVi = N(V S , V S ), (52) 

order el : C V 2 + CtV s = 7V(F s ,Vi). (53) 

Equation (fBTJl is satisfied identically, as it should, because 
V S is the sum of two terms with wavenumber k and — k 
which are both zero eigenmodes of C. The components 
of the Vo(£) follow from the numerical results for <f> of 
section [n] We can now solve equation ((53|) to find V\ . 
This gives us 

V x = B 2 (T,t)e 2i( - kz+UJrt '>- 2£t V 1 (0 + C(T)Vo + c.c. (54) 

where we have to find the vector V\(y) numerically, be- 
cause N(V) contains only quadratic terms. The second 
term is allowed with arbitrary C, since Vq is a zero mode 
of the linear operator Co ; it will not be needed in the sub- 
sequent analysis. In principle, there would also have to 
be an additional k = term, but this term is identically 
zero for the case of the asymmetric boundary conditions 
of interest here. Once V\ is known, we can proceed to 
Eq. ((53(1 : in this equation, the operator Ct works on 
B(T) and its complex conjugate only: we can write the 
equation as 

C V 2 + (Vo C T B + V * C T B*) = N(V , Vx). (55) 

This equation determines the time derivative of B 
through a solvability condition: since the operator Co 
has a right zero mode, it can be solved if and only if the 
other two terms in the equation are orthogonal to the 
left zero mode of Co H3, El ■ This requirement gives us 
the desired equation of dJ5/dT. To make this explicit, 
we first have to define the adjoint problem and an inner 
product. 

We define a space of 5-dimensional vector functions of 
the three variables £, z and t, 

«-> I/:/." •/.•'•:/ f(t, t) e *C— -')}. (56) 
The components of the functions / satisfy the physical 
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boundary conditions discussed in section ^ Further- 
more, an inner product / on f2 is defined: 



1 



2tt 



2 7-i « Jo w 



2tt 



r Jo 



(57) 



1=1 



On this space of functions, an adjoint operator £q is de- 



fined such that 



I(w,C Q f)=I{Clw,f) 



(58) 



for every function / in J7. Because the adjoint operator 
is obtained through partial integrations with respect to 
£, the requirement that does not yield any boundary 
terms from these partial integrations yields the appropri- 
ate boundary conditions for the functions g in the adjoint 
space. We will state these explicitly for our case in sec- 
tion nnu 

Let Wq be the zero mode of Cq. The solvability condi- 
tion applied to 155|) then becomes 

I(W , N(V Q S , Vi) - VoC T B - V *B*) = 0. (59) 

Let us focus on the term with B] if we write 

Wo = 9m(y,T,t)e lm ^ kz ^ (60) 

mgZ 

we see that we only need the term with m = 1, since 
the inner product in (|59|l is seen to vanish for all other 
terms after the z-integration is performed. For the same 
reason, we only need the terms which are proportional 
to e^z-^rt) f rom tv(Vo,Vi). Using this we obtain the 
following equation for the derivative of B: 

1 r 1 

- J d£ (W$, -C T B + \B\ 2 BN(V \ Fi)) = 0. (61) 

The complex conjugate of this equation is obtained from 
analyzing the term with B* in (|59|l . 

In (|59|) , we have used the approximation that e is suf- 
ficiently small that we are allowed to write 



r 



dt B(T)e~ £t (•) = B(T)e~ 



dt (■). 



(62) 



This is nothing but the usual adiabaticity assumption. 
As we have shown in section ^ this approximation is 
justified for large Wi. We can summarize the Eq. (|61(l in 
the form 



d T B = c \B\ 2 B, 



(63) 



where c 3 is a complex quantity which is obtained from 
working out the two inner product terms. We can now 
translate this result back into the lowest order nonlinear 



equation for the amplitude A. Upon combining (|49() , and 
(|63[) . we finally get 



d t A 



-eA + c 3 \A\ 2 A, 



(64) 



which is nothing but Eq. JBJ of the introduction in di- 
mensionless units. As discussed there, the sign of the real 
part of C3 determines whether or not we are dealing with 
a subcritical bifurcation. 

In the following sections the boundary conditions of 
the adjoint problem are derived, and we then proceed to 
solve adjoint problem and to determine c 3 . 

The framework laid out above can be extended rather 
easily to the case of the cylindrical tube. In that case we 
consider axially symmetric perturbations only, so that 
both Tg r and Tg z vanish; note, however, that rgg is non- 
vanishing in this case. The vector V therefore now has 
six components, Sv r , Sv z , Sr rr , 5T rz , Srgg and 6t zz . Apart 
from this change, the structure of the expansion is es- 
sentially the same, except for trivial changes like the fact 
that the first integration in l|57|l should now be taken over 
the two-dimensional scaled radial coordinate £. 



B. The adjoint operator and associated boundary 
conditions 



In this section we will calculate the components of the 
adjoint operator C\ using the defining equation l|58|) . We 
will follow again the planar case, and indicate the ma- 
jor changes for the cylindrical case at the end an in 
the appendices. Writing V — (vx,...,vs) and W = 
(wi, . . . , W5) we have 



C V 



d v vi + d z V2 



d y d z v 3 + (d 2 - 

CV! 

Dvi + Ev2 
Fvi + Gvi - 



d 2 )v4 - d y d z v 5 



) 2 
-Av 3 
- Bv 3 - 

2Bv 4 



Av A 
-Av 5 



(65) 



The various functions and coefficients in this expression 
are given in Appendix 1X1 We will illustrate the structure 
of the calculation by analyzing two terms of I(W, CoV), 
and simply state the results for the other terms. One 
term we get is w 3 Av 3 : 

jw* 3 Av 3 = J w *(l + d t + ^v°(0d z }v 3 , 



l-d t --v° z (0d z )w* 3 v 3 



+ 2 J9 z (w* 3 v 3 ) + Jd t {wlv 3 ). (66) 

Both integrals on the last line vanish, the first one be- 
cause the integrand is a partial z-derivative of a term 
which is periodic in z, the second one because the inte- 
grand is a partial i-derivative of a term which is periodic 
in t. Thus we simply obtain 



w 3 Av 3 = 



z w 3 v 3 . 



(67) 



13 



In short: we pick up a minus sign for every partial in- 
tegration and in performing the partial integrations we 
get, in general, boundary terms which have to vanish. 
In the above example, the boundary terms trivially van- 
ished because of the periodicity of the terms with respect 
to t and z, but the partial integrations with respect to 
£ do not automatically vanish. In particular, from the 
various terms we get the following boundary conditions: 

w* 3 vi(l) - = 0, 

V)lvi(l) + W%V!(-1) = 0, 

wlv 2 (l) ~ wlv 2 (-l) = 0, 
w* 5 v 2 (l) + w*v 2 (-l) = 0, 

-<ui(-l) = 0. (68) 

and 

w%d x v 3 (l) -w%d x v 3 (-l) = 0, 
w 2 d y V4(l) — w 2 d y V4(— 1) = 0, 
dyW^v^l) — dyW^v^— 1) = 0, 
w* 2 d z v 5 (l) - w* 2 8 z v 5 (-l) = 0. (69) 

The first set of conditions l|68|l is always satisfied, because 
V satisfies the boundary conditions of the original prob- 
lem: vi(±l) = v y (±l) = and u 2 (±l) = v z {±l) = 0. 
By setting w 2 = at £ = ±1, we have only one condition 
left of the second set: 

C (1K« 4 (1) - d £ w 2 (-l)* Vi (-l) = 0. (70) 

In order to understand how many boundary conditions 
we need to impose, it is good to realize that the operator 
C actually works on physical, and hence real functions; 
thus, we need the real part of this combination to be zero. 
Thus, we can still choose the phases of both amplitudes: 
we use this freedom to set the phase difference between 
the adjoint solution Wq and the solution V such that the 
combination to^ 174(1) becomes purely imaginary. This 
shows us that if we impose 

Re{d ( w 2 \ lVi (l) - ^m_ii7 4 (-l)) = (71) 

all conditions are satisfied for physical functions. We can 
give all wt in terms of w 2 and its derivatives. This gives 
a fourth order equation for w 2 which we solve using a 
shooting method. We will discuss the appropriate sym- 
metry conditions in the next section. 

In the cylindrical case, we need w 2 to be bounded. In 
the present case, this can be achieved by putting u;2(0) = 
0, w 2 (0) = 0. These conditions are satisfied if we expand 
w 2 (£) = ai£ + |a 3 £ 3 , for < x < 0.01. Furthermore 
we can set a± = 1 to eliminate the degree of freedom we 
have. This leaves us with one shooting parameter, 03. 
We again have the boundary condition that w 2 vanishes 
on the wall, ^2(1) = 0. We choose 0,3 such that this 
condition is satisfied. 



C. Symmetry of the solution and the adjoint 
solution 

As we discussed in section section [n] the linear planar 
mode that we investigate is asymmetric, which means 
that its symmetry is 

V(£) = (odd, even, even, odd, even) (72) 

which can be verified, see Appendix [BJ For the adjoint 
solution we have also have two zero modes with different 
symmetry, 

Wq 1 (£) = (even, odd, even, odd, even), (73) 

and 

Wq (£) = (odd, even, odd, even, odd). (74) 

It turns out that for our choice of the linear mode 172fl . 

(2) 

the second mode W does give a solvability condition. 
This can be seen as follows. Only the last three com- 
ponents of the vector V contain a derivative dr- Thus 
only the overlap of the third, fourth and fifth component 
of Wq and Vq come in; because of the symmetry of the 
opposite symmetry of the components, all these integral 
vanish. The same holds for the iV-term, and hence the 
solvability condition is trivially satisfied. As a result, 
Wq 1 ^ is the adjoint zero mode which gives the nontrivial 
boundary condition. The oddness of the component w 2 
is then guaranteed by taking 

w 2 (0) = 0, w' 2 (0) = l, w' 2 '(0) = 0. (75) 

These conditions, together with the boundary condition 
at £ = ±1 completely fixes the adjoint zero mode. As 
with the linear problem discussed in section [HJ we solve 
the differential equation together with the boundary con- 
ditions with a shooting method. 

In the cylindrical case, we have just one single mode, 
also for the adjoint problem; the boundary conditions 
that we already derived above for the cylindrical case are 
completely analogous to the above boundary conditions 

D. Numerical Results 

In summary, the nonlinear term C3 whose real part gov- 
erns the weak nonlinear stability is obtained numerically 
as follows. We first solve the fourth order equation for the 
stream function (f>o, derived in section [H] This gives us 
the components of the vector Vb(£). A second routine of 
this program generates the term N(Vo, Vo) too. A second 
program calculates the stream function (f>\ of the inhomo- 
geneous equation 15211 from which we obtain the compo- 
nents of the vector V\(y). A third program solves the ad- 
joint problem and gives the vector Wq. Because we have 
Wq ~ e i(kz-u r t) y ^ e i(kz—u r t) an( j y 1 ^ £ 2i(kz+uj r t) 
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FIG. 10: The value of the critical value A, as a function of FIG - 11: The value °1£ corresponding to the minimum value 
the wavenumber k for three different values of the Weissenberg of the curves m Fl S- HH1 as a function of Wi , and the values of 
number Wi for the case of the cylinder. fc at wnicn tne amplitude A is 1.1 times the minimum value. 

The jumps of the curves are due to fact that a curve of A c as 
a function of k has several minima. 



we have trivial z and t integrations in the inner product 
(|61[) . The ^-integration is done numerically from to 1 
because of symmetry. This then gives us the coefficient 

C3- 

We first discuss the cylindrical case. In Fig. ^| we 
plot the value of the critical amplitude A c (k), determined 
from C3 by Eq. 0, as a function of the wavenumber k for 
three different values of the Weissenberg number. The 
curves illustrate that for Wi > Wi c there is a band of 
wavenumbers where Rec3 > and that hence there is, in 
our approximation, a critical amplitude beyond which the 
flow is unstable. The sharp rise of the critical amplitude 
at the edges of the bands shown in Fig. ^| are close to 
the edges of the band where Rec3 — > 0. Note also that 
for Wi = 7.5, the critical value still has a rather sharp 
minimum near k = 4, but that for increasing values of Wi 
the bottom of the band flattens rapidly. For decreasing 
values of Wi, especially when Wi approaches about 5, the 
band likewise sharpens; we will see this from a different 
perspective below. 

Fig. ^]] also shows the behavior of the critical ampli- 
tude A c (k) as a function of k shows rather complicated 
structure. For Wi = 7.5 there is a plateau in the critical 
value of A c {k) around k — 2.7; upon increasing Wi, this 
plateau shifts and disappears, whereas the absolute min- 
imum of the curve shifts to smaller k values while a new 
minimum develops at larger k. Already at Wi = 8.25, 
the two minimum almost correspond to the same values 
of A Cl but for slightly larger Wi the minimum at larger 
k value becomes the absolute minimum. This is further 
illustrated in Fig. ^2 where we plot the values of k cor- 
responding to the absolute minima of the curves, as well 



as the values corresponding to an amplitude 1.1 times 
the minimum value, as a function of Wi. The figure illus- 
trates that between Wi ss 6.5 and Wi « 8, the minimum 
of the curve shifts to smaller k values, but that around 
Wi « 8 a local minimum at higher fc-values becomes the 
absolute minimum. 

As explained in section lll B 41 our normalization is such 
that A is the maximum perturbation in the shear rate 
at the wall over one wavelength, divided by the shear 
rate at the wall of the unperturbed parabolic profile [Sec 
Eq. ©]. Upon increasing Wi, the minimum values of 
the curves, which were already plotted in Fig. |3 quickly 
decrease. As we already mentioned in the introduction, 
we have also analyzed the critical value for the relative 
shear stress perturbation at the wall beyond which the 
flow is unstable [see Eq. ©]. The data for this ratio as 
a function of k for the same values of the Weissenberg 
number as in Fig. 1101 are shown in Fig. 1121 There are 
several important things to note about this figure. First 
of all, the curves of the critical shear stress perturbation 
have just one minimum, contrary to those for the criti- 
cal shear perturbation, and secondly, the values for the 
critical shear stress perturbation are typically a factor 
ten smaller than those for the critical shear. As Fig. IT2l 
shows, for Wi = 9.5 a perturbation of about 1% in the 
wall shear stress is sufficient to render the flow unstable. 
This is why in Fig.[3]the values of the critical shear stress 
amplitude (the values corresponding to the minima of the 
curves in Fig. I12f> are multiplied by 10 to draw them on 
the same scale as A c . Finally, we note that the edge of 
the band of unstable shear stress perturbations r is the 



15 




FIG. 12: The values of the dimensionless shear stress r be- 
yond which the flow is unstable in the cylinder, as a function 
of k for the same three values of Wi as in Fig. 1101 



same as the edge of the band of unstable shear pertur- 
bations — this is simply due to the fact that the edge of 
the band is marked by the point where Rec3 = 0. 

In Fig. 1131 we plot the band in which modes are nonlin- 
early unstable beyond some critical value as a function of 
Wi. This figure clearly illustrates that the width of the 
band vanishes at Wi = Wi c ps 5, since below this values 
Re C3 < for all k. This figure is somewhat reminiscent 
of the so-called "Busse balloon" of Rayleigh-Benard con- 
vection |3fjl |. but one should keep in mind that the inter- 
pretation is slightly different as we are dealing with a sub- 
critical (inverted) bifurcation. Thus, the regions marked 
"Stable" indicate the regions in the diagram where the 
perturbations with a wavenumber k in that region are 
nonlinearly stable in our approximation. The basic Pois- 
seuille flow profile is, however, nonlinearly unstable at all 
Wi c to modes whose wavenumber lies in the band marked 
"Unstable". In reality, we expect that the flow pattern 
will settle to some kind of stable nonlinear behavior dom- 
inated by modes in this band. 

Although the band of unstable modes widens toward 
small fc-values for Wi between 6 and 7.5 — a feature 
which is related to the plateau in Fig. ^|for Wi = 7.5 — 
at flow rates (Weissenberg numbers) about 50% beyond 
the critical value, the band of unstable wavenumbers 
ranges from just over 2 to a value somewhat larger than 
4. Although with our expansion to cubic order in the am- 
plitude we can not probe the stable nonlinear patterned 
flow regime nor the nonlinear selection of the wavenum- 
ber (it need not correspond to the most unstable mode), 
it seems reasonable to assume that the wavenumbers of 
the pattern in the tube will lie in the range identified by 

Fig. ma 



FIG. 13: Plot of the width of the fc-band where the corre- 
sponding modes render the basic cylindrical Poiseuille flow 
nonlinearly unstable. The dashed line indicates the fc-value 
of the mode which is most unstable to shear perturbations, 
i.e. the minimum of the curves plotted in Fig. 1101 

Because near Wi c the critical amplitude is relatively 
large, and again because we have only expanded up to 
cubic order, one has to interpret our results near Wi c 
with caution. Nevertheless, our results appear to identify 
the location of the saddle-node bifurcation to nontrivial 
flow patterns near Wi c ps 5. Since the band of unstable 
modes in this limit is very small, it appears likely that if 
one follows the nontrivial flow pattern down to this value, 
its wavenumber is should be close to the value where the 
band in Fig, lipidoses. 

Wi ps Wi c ps 5 : => k ps 3.75, (76) 

since the basic profile is then nonlinearly stable to pertur- 
bations with wavenumber that differs significantly from 
this value. In dimensional units, this implies for the wave- 
length A of the pattern 

Wi Pa Wi c ~ 5 : A ps 27T.R/3.75 m 1.7 R. (77) 

It is important to keep in mind that k is the wavenum- 
ber inside the cylindrical tube. If the (near)periodic flow 
and stress patterns inside the tube indeed cause the sur- 
face undulations that are the first stage leading to melt 
fracture outside the tube at larger flow rates, then one 
should keep in mind that k is not the wavenumber of 
these surface undulations. For, the polymer extrudate 
swells upon leaving the tube, and the flow velocity profile 
of the polymer becomes essentially constant after leaving 
the tube. A better way to compare our theory with mea- 
surements on the extrudate is therefore to compare the 
(dominant) frequency of the undulations: since no os- 
cillations will disappear at the outlet of the tube, the 
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frequency measured inside the tube must be equal to the 
frequency with which the the extrudate width oscillates 
after flowing out of the tube ^j- Of course, the fre- 
quency of the oscillations on the nonlinear flow branch 
corresponding to the solid line in Fig. ^ can not be ob- 
tained precisely from our expansion; assuming that the 
nonlinearities do not change this frequency too much, we 
estimate it from our results for the linear modes. From 
our analysis of the linear eigenmodes in section[n]we find 
that the dimensionless frequency 1m to of the modes is to 
a very good approximation (better than to a percent or 
so) given by 



Imcj = fc-0.93/Wi. 



(78) 



This result implies that for large Weissenberg numbers, 
the phase velocity Imw/fc approaches 1. Since according 
to we measure velocities in units of v ma x, this means 
that the periodic stress pattern moves essentially with the 
maximum velocity of the flow in the large Weissenberg 
limit. This indicates that the linear mode more and more 
concentrated near the axis of the tube for large Wi . The 
above result gives as an estimate for the frequency / in 
dimensional units 



/ = v max (k - 0.93/Wi)/(27ri?) 



(79) 



Finally, using that we expect the transition to occur at 
Wi c « 5 and k k, 3.75, and that for the unperturbed 
(parabolic) profile of a UCM fluid w max = 2v av with v av 
the average flow velocity in the tube, our estimate for the 
frequency at the transition becomes 



at transition: / w 1.13i> av /-R. 



(80) 



ft is important to realize that this estimate is based on 
the assumption that the frequency (or phase velocity) 
does not get renormalized significantly by the nonlinear- 
ities. This is maybe not unreasonable near threshold, 
where the amplitude of the periodic modulation of the 
pattern will be smallest (keep in mind that since we are 
dealing with a subcritical bifurcation, the amplitude near 
threshold remains finite). 

Especially further above threshold, we expect the 
renormalization of the frequency to be significant. Al- 
though this can not be calculated from our first non- 
trivial nonlinear term in the amplitude expansion, we 
can reasonably estimate the frequency at values of Wi of 
the order of 50 to 100% above threshold, say, as follows. 
As we noted above, in this range, the band of unstable 
wavenumbers k ranges from just over 2 to about 4.5 (See 
Fig. ll4fl . It is unlikely that the range of wavenumbers will 
change drastically due to nonlinear interactions, since in 
according to our calculations outside this range both the 
linear and the cubic term in the amplitude expansion are 
stabilizing. On the other hand, the frequency of oscilla- 
tions, if one studies the pattern in a fixed lab frame, will 
get strongly renormalized: once the pattern gets well de- 
veloped, we expect (but have no proof) that the nonlin- 
early modulated profile moves with the average speed v av 
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FIG. 14: Plot of the critical amplitude A c as a function of k 
in the planar case for three different values of Wi. 



rather than with the maximum speed of the unperturbed 
parabolic profile (keep in mind that in linear order, the 
periodic linear eigenmode does not affect the average 
speed, as it averages out to zero over one wavelength, but 
that this does not remain true in nonlinear order). Which 
particular wavenumber will be selected nonlincarly (or 
whether in fact a well-defined sharp wavenumber will be 
selected nonlinearly in a carefully controlled experiment) 
we do not know, but with the assumption that it lies in 
the unstable band fc m ; n < k < fc max and that it moves 
with the average speed, we get the following estimate for 
the dimensional frequency |33J 



2ttR 



</< 



27ri? 



which gives 



0.3 v a 
R 



<f< 



0-72 u a 
R 



(81) 



(82) 



We finally turn to a brief discussion of our results for 
the case of a planar slit. The main difference between the 
cylindrical and planar case is that in the latter case, the 
coefficient Re C3 is always found to be positive. This is 
also illustrated by Figs.^Jand^l which show that in the 
case of the planar slit there is no finite band of nonlinearly 
unstable wavenumbers for the critical shear rate (Fig. ll4|) 
and critical shear stress fFig. H5[l as a function of k. Note 
also that in both cases these curves just have a single 
minimum, contrary to what we found for the cylindrical 
tube. As a result of this, the wavenumber corresponding 
the minimum of the amplitude shear rate critical value 
A c now decreases smoothly with Wi, see Fig. 1161 

The absence of a finite band of wavenumbers and of 
a value of Wi below which Re C3 < for all k (as in the 
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FIG. 15: Plot of the critical amplitude r c as a function of k 
in the planar case for three different values of Wi. 



FIG. 16: The value k c of the wavenumber corresponding to 
the minimum of A c as a function of k for the case of planar 
Poiseuille flow. 



cylindrical case), appears to have fewer practical implica- 
tions than one might expect at first sight. After all, the 
overall behavior of the minimal value of the critical am- 
plitudes as a function of Wi shows quite the same behav- 
ior as in the case of the cylinder, compare Figs. El and 
upon decreasing Wi, the critical values rapidly decrease 
below Wi ps 5. As we have stressed before, our expansion 
ceases to be valid in this regime where A c becomes of 
order unity. We expect that in reality there is also a true 
saddle-node bifurcation point where the branch ends in 
the case of planar Poiseuille flow; possibly, our results in- 
dicate that in the planar case the corresponding critical 
value Wi c is lower than in the cylindrical case. However, 
our results do imply that it is difficult to make a pre- 
diction for the wavenumber near this (presumed) critical 
value than in the case of the cylinder. Nevertheless, since 
the wavenumber of the mode whose threshold to nonlin- 
ear instability is smallest is about 1.8, it seems reasonable 
to expect for the dimensional wavelength near threshold 



Wi re Wi, 



A pa 27rd/1.8 « 3.5 d 



(83) 



where the width of the slit equals 2d. Like for the case of 
flow in a pipe, according to our linear results the phase 
velocity is very close to the maximum velocity at the 
center line between the plates, so the above result for the 
most likely wavelength near threshold yields a frequency 
of about 2.7 v w /d. However, since we do not find a finite 
band of unstable modes above threshold, it is difficult 
to determine the range of frequencies expected further 
above threshold, even though we do expect here too that 
the patterns move roughly with the average velocity. 



IV. DISCUSSION AND OUTLOOK 

In this article we have shown that for the simplest 
model of a polymer fluid with normal stress effects, the 
UCM model, Poiseuille flow through a planar channel 
or cylindrical tube becomes weakly nonlinearly unstable 
for Weissenberg numbers somewhat larger than unity. 
Stated differently, since the UCM model only includes 
the essential normal stress effect, we find that the non- 
linear flow instability is characterized by the Weissenberg 
number only, and the phenomenon appears to be very ro- 
bust in that almost any more complicated polymer fluid 
model that includes normal stress effects will exhibit the 
same instability in the same ran ge o f Weissenberg num- 
bers. We presented evidence in |9Ul0j that this instability 
yields an intrinsic route to melt fracture behavior in the 
absence of other mechanisms such as stick-slip phenom- 
ena. 

One should also keep in mind that our expansion is 
only carried out to lowest order in the nonlinearity, so 
one may wonder about the robustness of these results 
as long as higher order terms in the expansion are un- 
known. Investigations of these issues for Couette flow 
and Poiseuille flow are underway and will be reported in 
due course. 

The critical Weissenberg number we find is a factor 
10 larger than the values for which Atalik and Keunings 
observed self-sustained quasi-periodic oscillations in 
their numerical simulations of Poiseuille flow. The two 
results are not necessarily inconsistent however: while 
all of our calculations apply to the UCM model, their 
simulations were done for the Oldroyd-B model, with a 
viscosity ratio of 10 -3 and Reynolds number Re = 0.1. 
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Moreover, Atalik and Keunings added a stress diffusion 
term to their equations to improve numerical stability. 
Further analysis is clearly necessary to investigate the 
effects of these differences — a detailed comparison of 
the simulations with the amplitude expansion results for 
one and the same model, is called for. 

Just above the onset of a well-defined linear (or "su- 
percritical" ) instability threshold, the wavelength of the 
pattern is well-defined: upon approach of the threshold, 
the wavenumber of the pattern approaches the wavenum- 
ber k c at which the instability sets in. Here, however, 
it is more difficult to draw sharp conclusions about the 
wavelength of the pattern close to onset. First of all, 
our expansion can not be fully trusted quantitatively for 
Weissenberg numbers of the order of Wi c , as the damping 
rate is large there, not small as is needed in our method. 
This technical caveat aside, one should keep in mind that 
once the instability sets in the nonlinear flow behavior 
which will develop can not be addressed by our expan- 
sion method. Hence it is difficult to rule out the possibil- 
ity of nonlinear interaction terms changing the velocity 
of the pattern in the die to a value different from the 
one we have assumed, namely the average flow velocity 
v a . If on the other hand the flow pattern stabilizes in 
some weakly nonlinear regime, then one would expect its 
wavenumber to be close to the onset value we calculate 
and the frequency to be close to the one we have esti- 
mated in Eq. (|80l) . As a rule of thumb, the wavelength of 
the undulations of the extrudate is typically about twice 
the diameter of the die. 

In conclusion, our nonlinear analysis establishes the 
nonlinear flow instability essentially beyond reasonable 
doubt and predicts onset values which are consistent with 
those reported experimentally 0, H, El El • Moreover, 
the hypothesis of highly similar subcritical behavior in 
Taylor- Couette geometries El EJ and Poiseuille flow 
(and hence possibly melt fracture) is fully confirmed by 
our calculations. Indeed, the similarity between the flow 
field perturbation shown in Fig. [S] and the roll- type pat- 
tern which Groisman and Steinberg have argued gives 
rise to the subcritical nature of the instability [2(| in 
Taylor- Couette flow is striking. From this perspective, 
the only difference between the two cases is that the gen- 
eral mechanism |0, El 03 that in visco-elastic flows the 
curvature of the streamlines makes the flow linearly un- 
stable is operative in Taylor- Couette cells but obviously 
not in Poiseuille flow. Nonlincarly, the two flows ap- 
pear to be much more closely connected. From this per- 
spective, the evidence for "turbulence without inertia" 
[34l l35| as a result of these elastic effects is intriguing. 
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APPENDIX A: EXPLICIT EXPRESSIONS FOR 
THE OPERATOR C AND THE NONLINEAR 
TERM TV 

In this section we give the expressions for the operators 
C and N from the equation 



CV = N(V, V) 
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in the planar case and the cylindrical case in dimensional 
form. 

For the planar case we have 
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where the operators A — G are defined in the following 
way: 
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For the cylindrical case we get 



and the vector N has components 
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where the operators A — M are defined in the following 
way: 
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APPENDIX B: EXPLICIT EXPRESSIONS FOR 
THE COEFFICIENT ft 



The coefficients of the equation for the streamfunction 
in the planar case: 
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where the complex dimensionless coefficient c = ui/k is 
used by Rothenberger et dl @. 

The coefficients of the streamfunction in the cylindrical 
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0o := fc 3 (fc-4iC(OS 3 C 2 +4fe[SYC(0 2 
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-3 + 4S 2 £ 4 C(£)fc 2 + 2i fc 3 £ 4 S 
-2C(0 2 S 2 C 4 fc 2 + 4* fc 3 ^ 6 C(0 2 S 3 
-4ifc 3 £ 6 C(£)S 3 )/(£ 3 ), (B7) 

02 := -(-3-4S 2 ^ 4 C(Ofc 2 + 2C(0 2 S 2 e 4 fc 2 

+2fc 2 £ 2 + 2S 2 £ 4 fc 2 )/(£ 2 ), (B8) 

3 := 2(-*fce 2 S-l + *Se 2 C(e)fc)/^ (B9) 
(7(0 := l/(l + i/2fcS(l-£ 2 -c)). (BIO) 
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